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Kozai-Lidov oscillations of Jupiter-mass planets, excited by comparable planetary or brown dwarf 
mass perturbers were recently shown in numerical experiments to be slowly modulated and to 
exhibit striking features, including extremely high eccentricities and the generation of retrograde 
orbits with respect to the perturber. Here we solve this problem analytically for the case of a test 
particle orbiting a host star and perturbed by a distant companion whose orbit is eccentric and 
highly inclined. We give analytic expressions for the conditions that produce retrograde orbits and 
high eccentricities. This mechanism likely operates in various systems thought to involve Kozai- 
Lidov oscillations such as tight binaries, mergers of compact objects, irregular moons of planets 
and many others. In particular, it could be responsible for exciting eccentricities and inclinations 
of exo-planetary orbits and be important for understanding the spin-orbit (mis)alignment of hot 
Jupiters. 



A Keplerian orbit weakly perturbed by a distant orbit- 
ing mass may exhibit long term, large amplitude cycles 
in which the eccentricity and inclination change periodi- 
cally [2. 3]. These so-called Kozai-Lidov cycles are owed 
to the orbit-averaged quadrupole potential of the per- 
turber. The high eccentricity excited from initially nearly 
circular orbits by this mechanism is suggested to play an 
important role in the formation and evolution of many 
astrophysical systems [e.g. I3H81. Hi!] . 

Kozai-Lidov oscillations of Jupiter-mass planets, ex- 
cited by comparable planetary or brown dwarf mass per- 
turbers were recently shown in numerical experiments to 
be slowly modulated and to exhibit striking features, in- 
cluding extremely high eccentricities and the generation 
of retrograde orbits with respect to the perturber [ll[ 
(see also [HI). Those authors attributed the slow mod- 
ulation of the Kozai-Lidov cycles to the relative prox- 
imity (including the octupole potential of the perturber) 
and comparable mass of the perturber to the perturbed 
planet, and the high eccentricities to chaotic evolution. 
These effects were argued to play a major role in affect- 
ing the properties of hot Jupiters (extra-solar, Jupiter 
mass planets with very short periods < 10 days) and in 
particular in explaining the large fraction of hot Jupiters 
recently found to have a retrograde orbit with respect to 
their host's spin. 

In this letter we show that the Kozai-Lidov cycles are 
slowly modulated non-chaotically and quite simply by 
the octupole potential of the perturber, which is non- 
vanishing when the perturber's orbit is eccentric. This 
slow modulation can excite extremely high eccentricity 
and inclination of an initially nearly circular Keplerian 
orbit. A particular consequence is that the orbit plane 



"John Bahcall Fellow, Einstein Fellow 
tSagan Fellow 



can flip to retrograde with respect to the total angular 
momentum of the system These effects occur in 

the test particle approximation in which the mass of the 
perturber is much bigger than that of the planet, and will 
thus be important for stellar binary perturbers of planets. 
We describe this long term evolution of the Kozai-Lidov 
cycles by deriving and analytically solving the effective 
equations, averaged over the Kozai cycles. 

Secular Equations Consider a test particle on a Kep- 
lerian orbit (semi-major axis a and eccentricity e) subject 
to perturbation by a distant mass Af por on an orbit (a peT , 
e pcr ) around the same central mass M. The coordinate 
system is defined using the perturber's orbit, with the 
z-axis chosen to be in the direction of the angular mo- 
mentum vector and the x-axis pointing to the pericenter. 
It is useful to parametrize the test particle's orbit by two 
dimensionless vectors: j = J/yGMa, where J is the spe- 
cific angular momentum vector and G is the universal 
constant of gravitation; e, a vector pointing in the direc- 
tion of the pericenter with magnitude e. The orientation 
of j is defined by the inclination i with respect to z and by 
the longitude of ascending node f2 (angle between z x j 
and x), as follows j = j( sin i sin f2, — sinzcosJl, cosi). 
Usually, the orientation of e is set by additionally speci- 
fying the argument of pericenter u> (angle between e and 
z x j). Here we define the orientation of e by the co- 
latitudinal angle < i e < ir (angle between z and e), 
and longitude Q e (angle between the projection of e on 
the xy plane and x) by 

e = e(sinz e cosJ7 e5 sini e sinf2 e , cosi e ). (1) 

It turns out that for the cases considered, £l e is slowly 
varying and is useful for describing the long term behav- 
ior of the system. 

The secular orbital evolution of the test particle is de- 
termined by double time-averaging the perturbing po- 
tential $p Cr over the orbital periods of the test parti- 
cle and the perturber. The averaged potential expanded 



2 



to the octupole order (3rd order in a/a peT ) is given by 
(3> P er) = = *o(^Quad + eoct^Oct), where the dimen- 



sionless averaged potential 4> is expressed as the sum of 
two components (quadrupole and octupole), 
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and the normalization parameters are 
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In the secular approximation, a and (f> are constant 
with time while j and e evolve according to the following 
equations of motion [t| [T2, [H| , 



dr 



= jxVj 
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= jxV e ^ + exVj^, (6) 



where t = t/t scc and t scc = V GMa/$>o is the secular 
timescale. Physical solutions are restricted to those sat- 
isfying the physical constraints j 2 = |j| 2 = 1 — e 2 , and 
e ■ j = 0. 

Kozai-Lidov Cycles When expanded to the 
quadrupole order only (i.e., eoct = 0), the averaged 
perturbing potential is axisymmetric. As a consequence, 
j z is conserved and the equations of motion are in- 
variant under rotational transformations around the 
z axis. In this case, e, i, uj and i e undergo periodic 
oscillations (Kozai-Lidov cycles), which are determined 
by the two constants of motion j z and cf> = 0Q ua d [e.g. 
U 0]. It is convenient to use the constant of motion 
C K = 
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\iz + \ -> which is given by 
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When Ck < 0, u> librates around 7r/2 or — ir/2 (Kozai- 
Lidov librations), while for Ck > 0, uj varies monoton- 
ically with time taking all values from to 2tt (Kozai 
Lidov rotations). 

For rotations, e reaches minimum at uj = or 7r, im- 
plying 
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Ck (only when Ck > 0). 



(8) 



For librating solutions, e m - m is obtained at uj = ±7r/2. 
For any Kozai-Lidov cycle, maximum e is obtained at 
uj = ±7r/2, leading to 
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To fully specify the trajectory, the azimuthal position 
must be specified. We use the azimuthal angle il e , which, 
as shown below changes slowly for the regime j 2 < 1 of 
interest here. In fact, the equation of motion for tt e reads 



where fn = 3(8 



^ e = jzfn 
6/ sin 2 i e )/8. 
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Kozai-Lidov Cycles with j z = In this case, Eqs. 
©, and © imply that the torque is given by j = 
— (15/4)e z e x z, and is directed (up to sign) in the di- 
rection of j. This means that j moves on a straight line 
through the origin j = (at which e = 1) in the xy plane 
periodically with 
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Each cycle, when j crosses the origin, jumps by n (Q 
is constant otherwise). 

e is confined to the plane perpendicular to the line 
on which j moves and thus fi e is constant as long as e 
does not cross z. For e to cross z, e must equal e z so 
that Ck < 0. Hence, for rotating cycles (Ck > 0), e 
never crosses the z axis and f2 e is constant throughout 
the cycle. For librating cycles, e crosses z each cycle 
changing f2 e by ir. Note that the derivative of f2 e , given 
by Eq. (fTU)) , is zero for j z = except for a divergence at 
the point i e — 0, where e points in the z direction. 

Equations of motion We next study the evolution of 
the system including the contribution of the octupole 
part of the perturbing potential. Since the octupole con- 
tribution is small, we can assume that at any time, on 
short time scales (of order t sec ), the solutions are nearly 
Kozai-Lidov cycles. The properties of these cycles at 
any given time are determined by the values of j z and 
Ck at that time. Moreover, given that the total poten- 
tial (j) is conserved, we have to a good approximation 
<?^Quad = const, implying that 
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Thus, the problem can be reduced to finding j z (t)- 

The time derivative of j z arises from the octupole po- 
tential alone; using Eqs. ([2]), and (|6|), it is given by 
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We focus on Kozai-Lidov cycles with j 2 <C 1 and study 
the long term behavior of the system. Taking the lowest 
order terms in j z , we have 



jz = -e 0ct siii(S} e )f j (e,i e ) 



(14) 



where f 3 = (75/64)eshn e [| - e 2 (| - 7cos 2 i e )] . 

For rotating cycles, fi e varies slowly implying that j z 
changes slowly. For librating cycles, f2 e changes by nearly 
7r each cycle, implying that j z goes to approximately —j z 
and the contribution to j z vanishes to zeroth order in j z 
(see example in figure [3]). Below we focus on rotating 
cycles for which j z can monotonically change over several 
secular timescales. 

Using Eq. (I12[) . and assuming that the initial condi- 
tions jz,a,CK,o are in the rotating zone (Ck,o > 0), the 
condition for rotation is — j Zima x < jz < Jz,max where 



= J2C Kt0 +j 2 



(15) 
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If j z crosses this border, j z changes sign during each cy- 
cle. For the examples of such cases we checked, after a 
few cycles j z moved away from this limit, back into the 
rotation region (see example in top panel of figureH]). We 
note that for some initial conditions, we found significant 
modulation (including flips) occurring for librating cycles 
on very long time scales t ~ eQ 2 t t scc , much greater than 
those studied here t ~ eQ* t t sec . 

Averaged equations We next average the equations of 
motion over the Kozai-Lidov cycle to obtain approximate 
equations that describe the long term behavior of the 
system. To lowest order in j z , eoct, we can average Eqs. 
(fl4t and (fTO]) over a cycle by taking the limit j z = in 
which tt e is constant and neglecting the deviation of fl e 
due to the octupole. The latter is important only during 
short episodes when \j z \ < eoct, in which j z can change 
considerably during one cycle and which are not resolved 
in the long term approximation. We obtain 



= jz (fa) 
jz = -eoct (fj) sin(fie), 



(16) 



where fi (i = Cl,j) are averaged over a Kozai cycle with 

Jz=0, 



(fi) = — l dtfi = — 

T~Koz Jj z =0 T Koz Jo 
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where tk oz = § dt is Kozai-Lidov cycle period at j z = 0. 
To evaluate the integral, note that e 2 (l — (5/2) cos 2 i e ) = 
Ck, which together with e 2 = 1 — j 2 and e^ nin = Ck 
allows a straightforward integration using (|lip . yielding 
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where K(m) = J*' (1 - TOsin 2 (6»))- 1 / 2 d6' and E(m) = 

fo^ 2 (l — m sin 2 (6)) 1 / 2 d9 are the complete elliptic func- 
tions of the first and second kind respectively. Note that 
(fi) are functions of Ck only since the Kozai-Lidov cycle 
over which the averaging is made has j z = 0. The upper 
panel of figure Q] shows plots of (/n) and (fj). 

Eqs. (TlTjjl . flTgl) and (Ti"2"j) form a closed set of equations 
for the slowly varying j z , Ck and tt e . An example of a 
numerical integration of these equations, compared with 
the results of a direct integration of the secular equations 
(|6]) is shown in figure [2] for eoct =0.01. As can be seen the 
approximate equations describe the long term evolution 
to a good approximation. 

These equations break down if \j z \ crosses the thresh- 
old Eq. (TT5j) . in which case fi e receives kicks, j z changes 
sign, j z moves to the rotation region after a few secu- 
lar time scales and the averaged equations become valid 





FIG. 1: Upper panel: (fj) (solid) and (fn) (dashed) from 
Eqs. (Tl8)) . Lower panel: F from Eq. ((20)) (solid line) and its 
quadratic approximation Eq. (|20[) (dashed line) 
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FIG. 2: Results of numerical integrations for initial conditions 
eoct = 0.01, uio = 0, fi = 7T,i = 80°, e = 0.1. The blue 
solid lines are the result of the direct integration of Eqs. ([6]) 
while the red dashed lines are the results of integrating Eqs. 
(|16p , (ll8[) and (|12[1 and using the pure Kozai relations to ex- 
tract e m i n , £max and i m in, imax- The two green horizontal lines 
in the top panel represent ±j max , given by Eq. (|15l) 



again. An example of such behavior is seen in fig 2) In 
this example, the effective equations Eqs. (fTB f . (TT5)) and 
(fT2)) were integrated only in the intervals were Ck > 
(dashed lines). 

Analytic solution Eqs. (|16p . (|T8l) and (TT21) are inte- 
grable. These equations have a constant of motion of the 
form 



C = F(C K )-e oct cos(n e ). 



(19) 
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FIG. 3: Results of numerical integrations for initial conditions 
identical to the parameters in figure [2] except for ujq = n/2 
replacing loq — 0, for which the Kozai cycles are rotating. The 
blue solid line is the result of the direct integration of Eqs. 
©. The two green horizontal lines represent ±j ma x, given by 
Eq.Jinj 
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FIG. 4: Results of numerical integrations for initial condi- 
tions eoct = 0.01, cjo = 0, Q.0 = 0,io = 88°, eo = 0.1. The 
blue solid line is the result of the direct integration of Eqs. 
([6]) while the red dashed line is the result of integrating Eqs. 
Q16 [I .(|18 P and (|12|l in the regions were Ck > (corresponding 
to \ jz | < Jmax, the region within the two green horizontal lines 
given by Eq.fl5|)) 



Indeed, using (fTp)) and (fT2"j) we find, C = 
e oct j z smtt e (F' (fj) - {fa)), and by setting 

P(r s f° K (/n) (c) , 

f{Ck)= L jim dc 

= 32^ f 1 - ~ 2 f *) dx (20) 
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we obtain (7 = 0. The numerical value of F{e^ nin ) as a 
function of e m - ln is shown in the bottom panel of figure Q] 
and tabulated in the attached file e2F.txt. 

Note that F diverges at e minj i n f = (4/11) 5 where 
{fj) = 0- F° r e min > emin,inf) the integration limits in 
Eq. ([20]) should be chosen differently. Here we focus 
on e m i n < e m ; ni i n f. F has a maximum F max w 0.0475 
at e min,m = ( 3 ~ 3a: ro )/(3 + 2x m ) « 0.112 where x m is 
the solution to the equation K(x m ) = 2E{x m ) (at which 
{hi) — 0). Near maximum, F can be well approximated 
by a quadratic expression, 

-^( e min) ~ ^iiiax — 2.67(e min — e minm ) , (21) 

(shown as dashed line in the bottom panel of figure [l}. 

This constant of motion holds most of the information 
about the system. 

Flip criterion We next use the constant of motion 
Eq. (fT9|) to derive a criterion of the initial conditions 
which is necessary in order that j z "flips" i.e. changes 
sign (implying i goes above 90° and the orbit becomes 
retrograde relative to the perturber). 

During a flip, j z = and Eq. (fl2|) implies that 
Ck = Ca'.o + 0.5j% . Given the constant of motion Eq. 
(JTnj), and that the term eoct cos fl e can change by at most 
2eoct, a required condition for a flip is that eoct > eoct.c 
where 

eoct, c = -max(|AF(x)|) (22) 

where x is in the range Ck,q < x < Ck.o + 0.5j^ and 
AF{x) = F{x) — F{Ck,o)- Given that F has one maxi- 
mum at Cx.m — ^min m ' ^ ne condition can be separated 
as follows: If Ck,o + 0.5j| < C'k >m , or Ck.o ^ Ck^tti^ so 
that the initial and final Ck are on the same side of the lo- 
cation of C K , m , eoct.c = \F(C K ,o + 0.5jf ) - F(C K ,o)\- 
Otherwise, if the initial and final Ck are on differ- 
ent sides of CK, m , the value of eoct is the larger of 
eoct.c = \F(C K ,o + 0.5jf )0 ) - F{C K ,o)\ and e 0c t,c = 
F mm — F(Ck,o)- The presence of F max thus creates a 
discontinuity in the flip condition (see figures [5] and [6]) . 
For cases where initially eo -C 1 implying that Ck *C 1 
and j zfi = cosi , and for j 2 zfl < 2e% dn m (i > 61.7°), Eq. 
(f2"2")l reduces to 

tOct,c= ^(icos^o). (23) 

This analytic theoretical line is shown in solid blue 
in the upper panel of figure [5j For comparison, the re- 
sults of numerical integrations of Eqs. ([6]) over 10/eoct 
secular times for corresponding eoct, io with ujq = 0, 
eo = 0.001 and fio scanned over — 2-7T are shown in filled 
red (flipped) and empty blue (no flip) circles. The ana- 
lytical curve describes the flip condition to better than 
10% for i > 80 deg, better than 20% for i > 70 deg, 
and to a factor less than 2 for i > 50 deg. The devia- 
tion at low inclinations is not surprising, given that our 
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FIG. 5: Upper Panel: The solid black line is the threshold, Eq. 
(1221 . for a flip given that we start with eo ~ (this reduces 
to (|23[) for i > 61.7°). Red circles are results of numerical 
integrations (eo = 0.001, cjq = 0, fio = — 2n) that had 
a flip while blue circles are simulations that didn't. Lower 
Panel: The solid black line is the flip threshold, Eq. for 
£Oct = 0.01. Red circles are results of numerical integrations 
(u>Q = 0, fio — — 2tt) that had a flip while blue circles are 
simulations that didn't. 

formalism assumes j z <C 1. It is encouraging that the 
overall behaviour is captured quite well for j z up to 0.5. 

The curve representing Eq. for eo c t,c = 0.01 on 

the e m i nj o,?o plane is shown in the lower panel of figure 
[5] (solid black). For comparison, the results of numeri- 
cal integrations for corresponding e m i n! o, io are shown in 
filled red (flipped) and empty blue (no flip) circles. 

The result of Eq. ([^)) for different values of e m i n is 
shown in figure [6] 

Maximal e and General Relativity ( GR ) precession A 
rough estimate of the typical maximal eccentricity ex- 
pected during an episode when j z crosses zero, can be 
obtained as follows and was confirmed in several numer- 
ical test runs with various parameters. The assump- 
tion is that at maximal e we have j ~ \j z \. Since the 
maximal e will be obtained in some arbitrary phase in 
the Kozai cycle during which j z crosses zero, we expect 
1 - eL = J 2 ~ Jz ~ (0.5j>koz) 2 ~ £oct (where we as- 
sumed that J7 e is not tuned to or n so that (sin 2 fl e <~ 1). 
In fact we expect a roughly uniform distribution of j z or 
[1 — max(e)] 1 / 2 around this range. General relativistic 
corrections cause the pericenter to rotate and can be in- 
corporated into the equations of motion by adding a term 
0GR = £gr/j to the total normalized averaged poten- 
tial where e G R = 3GM 2 /M pcl a- 4 al CI . This effect 
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FIG. 6: The theoretical threshold Eq. p2jl, for a flip 
for different initial e m i n (eo at u = 0). The black solid 
lines show the result of Eqs. (122I) and (|20[) . for e m i n = 
0.001,0.1,0.2,0.3,0.4,0.5 (lines in the top left corner are or- 
dered low eoct,c,io to high eoct, c ,io for growing e m i n ). The 
dashed blue curves are the corresponding results of Eq. (|22|) . 
with the approximate expression for F, Eq. (|2ip . 



suppresses the maximal eccentricity in the Kozai-Lidov 
cycles when eQ^/j ~ 1 [6;]; using the estimate above for 
e m ax) GR becomes significant once €qr ~ eoct- If this 
occurs, e can change its direction considerably during 
one Kozai-Lidov cycle, f2 e may change by ~ 7r, causing 
j z to change sign avoiding a flip. We verified this rough 
criterion for suppression of flips numerically. 

Kozai-migrated Hot Jupiters by distant, stellar mass 
perturbers Recently, numerical simulations by [Til ] 
showed that Jupiter-mass planets, perturbed by compa- 
rable planetary or brown dwarf mass perturbers, undergo 
slowly modulated Kozai-Lidov cycles, and exhibit strik- 
ing features, including extremely high eccentricities and 
the generation of retrograde orbits with respect to the 
perturber. Those authors attributed the slow modula- 
tion of the Kozai-Lidov cycles to the relative proximity 
and comparable mass of the perturber to the perturbed 
planet, and the high eccentricities to chaotic evolution; 
they suggested that these effects would not occur in the 
case of distant stellar mass perturbers. 

Our analytical and numerical analysis shows that these 
effects have a simpler explanation: the octupole pertur- 
bation alone modulates the Kozai-Lidov cycles, and it 
does so non-chaotically to excite the extremely high ec- 
centricities and the retrograde inclinations. This effect 
occurs already in the test particle approximation in which 
the mass of the perturber is much bigger than that of the 
planet, and will thus be important for stellar mass per- 
turbers. In fact, similar evolution occurs in other cases 
where there are small deviations from axisymmetry [lil ]. 

Consequently the distribution of orbital parameters of 
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Kozai-migrated hot Jupiters due to a stellar perturber 
may be significantly affected by the dynamics of the oc- 
tupole perturbations described in this Letter. 

Consider for example, the system parameters stud- 
ied numerically in [6], a — 5AU, a pcl = 500AU, 
M = Mpe r = M© ~ lOOOMj, where the contribu- 
tion of the octupole term was neglected. The per- 
turber eccentricity in this study was set to for conve- 
nience (which implies zero octupole contribution), but 
in reality is expected to have a wide distribution of 
values. Such a system has an octupole coefficient, 
eoct = aa per e per/(l — e p C i) ~ 0.01, secular time scale of 
t soc oc M 1 ' 2 M~^a- 2 al cr = 1.8 x 10 6 yr, and GR preces- 
sion coefficient of €gr oc M 2 M~ c ^a~ 4 ap OI . = 4.7 x 10~ 5 . 

We performed a few test runs with these parameters 
(including GR precession but not tidal dissipation) and 
found that extremely high eccentricities 1 — e max ~ 10~ 4 
can be reached for inclinations i > 80° (in reality, e max 



would be limited by other physical effects) . For the runs 
we made, a flip was suppressed due to the GR preces- 
sion, but for equally likely parameters with slightly closer 
perturbers, a pcr < 300 AU the effect of GR precession 
becomes negligible and flips are also attainable in accor- 
dance with the criterion Eq. (|22|) . 

A numerical investigation of this problem is published 
simultaneously by a different group in [Til l . 
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